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Self-organization of neural circuitry is an appealing framework for understanding cortical 
development, yet its applicability remains unconfirmed. Models for the self-organization of 
neural circuits have been proposed, but experimentally testable predictions of these models 
have been less clear. The visual cortex contains a large number of topological point defects, 
called pinwheels, which are detectable in experiments and therefore in principle well suited 
for testing predictions of self-organization empirically. Here, we analytically calculate the den- 
sity of pinwheels predicted by a pattern formation model of visual cortical development. An 
important factor controlling the density of pinwheels in this model appears to be the pres- 
ence of non-local long-range interactions, a property which distinguishes cortical circuits from 
many nonliving systems in which self-organization has been studied. We show that in the limit 
where the range of these interactions is infinite, the average pinwheel density converges to vr. 
Moreover, an average pinwheel density close to this value is robustly selected even for inter- 
mediate interaction ranges, a regime arguably covering interaction-ranges in a wide range of 
different species. In conclusion, our paper provides the first direct theoretical demonstration 
and analysis of pinwheel density selection in models of cortical self-organization and suggests 
to quantitatively probe this type of prediction in future high-precision experiments. 

1 Introduction 

Neuronal circuits in the mammalian cerebral cortex are among the most complex systems in na- 
ture. The biological mechanisms that contribute to their formation in early brain development 
remain poorly understood. However, it is unlikely that the precise architecture of mature cortical 
circuits can be attributed to genetic prespecification, since the number of genes in the genome is 
insufficientUJ. Instead, dynamical self-organization presumably plays a major role in shaping the 
architecture of neuronal circuits in the cerebral cortex. Dynamical self-organization is most thor- 
oughly described in non-living physical systems driven outside of thermodynamic equilibrium by 
external forcing. Whereas the emergence of structure is externally driven, the structures formed 
are primarily determined through interactions within the system itself [2,^. Neural circuits en- 
compass various positive and negative feedback loops, and these could well form the basis for 
cortical self-organization. However, evidence supporting this presumption is derived from theo- 
retical considerations IHISHSl rather than empirical observation. Models for the self-organization of 
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neural circuits have been proposed, but experimentally testable predictions of these models have 
been lacking. 

The system of orientation columns in the visual cortex is a paradigmatic system for studying 
cortical development and the role of self -organization in this process. Most neurons in the visual 
cortex respond selectively to a particular orientation of an elongated visual stimulus. Whereas in 
columns perpendicular to the cortical surface, neurons prefer similar stimulus orientations, the 
preferred orientation varies mostly smoothly[7J and repetitively across the cortical surface giving 
rise to a complex two-dimensional pattern called the map of orientation preference (Fig. la, b). 
Throughout the cortical map, there are point-like orientation singularities ||8l|9l[10l called pin wheel 
centers[llj at which all stimulus orientations are represented in circular fashion. Numerous stud- 
ies are consistent with the hypothesis that orientation maps develop through activity-dependent 
self-organization. They form in dark reared animals il2|. under substantial manipulation of visual 
input llT3llT4| . and even in auditory cortex when rewired to be driven by visual inputsfTS]. More- 
over, an analogy between cortical development and pattern formation appears plausible. Like in 
other systems where pattern formation has been observed HZHS], the orientation map arises prob- 
ably from an initially non-selective state, it exhibits a typical periodicity and a spatial extension at 
least an order of magnitude larger than the basic periodicity length. 

The conditions under which orientation maps can arise through self-organization have been 
thoroughly investigated theoretically l|T6l llTl [TSlI . A recent and highly promising approach stress- 
ing the analogy to pattern forming systems showed that pinwheels can be stabilized by activity- 
dependent long-range interactions lfT9| . A phenomenological order parameter field model based 
on the Swift-Hohenberg equation f20l was proposed in which orientation maps arise from a super- 
critical bifurcation of Turing-type. In this class of models, the stabilizing nonlinearity includes only 
key features of visual cortical organization and is constraint by biologically plausible symmetry as- 
sumptions. The model exhibits multiple structurally distinct quasiperiodic attractors resembling 
orientation maps in the visual cortex. 

The qualitative similarity of solutions of this model to orientation maps in the visual cortex 
appears promising. However, it is unclear at present weather the model accounts for aspects of 
cortical organization also quantitatively and whether this resemblance is not just superficial. Com- 
paring directly orientation maps in the model and the visual cortex would be difficult, because of 
the large number of possible map structures. Instead, it will be virtually unavoidable to take a 
statistical approach to this questions. As discrete entities, pinwheels can be detected, counted and 
localized with high accuracy. Here, we calculate in the model the average density of pinwheels in 
the limit of infinite interaction range. We show that this density is representative for a large regime 
of intermediate interaction ranges covering the estimated ranges in various mammalian species. 
Therefore, the quantity pinwheel density appears particularly well suited for testing for signatures 
of long-range dominated self-organization in experiment. In the following, we briefly describe the 
system of orientation preferences in the visual cortex and present the long-range self-organization 
model for its activity-dependent development. 

1.1 Orientation preference maps 

In the visual cortex, as in most areas of the cerebral cortex, information is processed in a 2-dimensio- 
nal (2D) array of functional modules, called cortical columns pTl |22) . Individual columns are 
groups of neurons extending vertically throughout the entire cortical thickness that share many 
functional properties. Orientation columns in the visual cortex are composed of neurons preferen- 
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Figure 1: Patterns of orientation columns and long-range horizontal connections in the primary vi- 
sual cortex of tree shrew visualized using optical imaging of intrinsic signals (modified from ||23|). 
a. Activity patterns resulting from stimulation with vertically and obliquely oriented gratings, re- 
spectively. White bars depict the orientation of the visual stimulus. Activated columns are labeled 
dark gray. The used stimuli activate only columns in the primary visual area VI. The patterns 
thus end at the boundary between areas VI and V2. b. The pattern of orientation preferences cal- 
culated from such activity patterns. The orientation preferences of the columns are color coded as 
indicated by the bars. A part of the pattern of orientation preferences is shown at higher magnifi- 
cation. Two pinwheel centers of opposite topological charge are marked by arrows, c. Long-range 
horizontal connections extend over several millimeters parallel to the cortical surface (tree shrew, 
superimposed on the orientation preference map). White symbols indicate locations of cells that 
were filled by a tracer (biocytin); labeled axons are indicated by black symbols. 

tially responding to visual contours of a particular stimulus orientation [Zl. In a plane parallel to 
the cortical surface, neuronal selectivities vary systematically, so that columns of similar functional 
properties form highly organized 2D patterns, known as functional cortical maps (Fig. H^, b). 

Experimentally, the pattern of orientation preferences can be visualized using optical imaging 
methods [8, 9J. Optical imaging of intrinsic signals is based on the fact that the optical properties 
differ in active vs. less active parts of the cortex i24| . This is utilized to record patterns of activity 
from light reflectance. In a typical experiment, the activity patterns -Efc(x) produced by stimulation 
with a grating of orientation 6^ are recorded. Here x = (x, y) represents the location of a column in 
the cortex. Using the activity patterns Ek{^), a field of complex numbers 2:(x) can be constructed 
that completely describes the pattern of orientation columns: 




(1) 
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The pattern of orientation preferences i9(x) is then obtained from 2:(x) as follows: 



^(x) = 2 o-T^aiz) 



(2) 



Typical examples of such activity patterns Ej. (x) and the patterns of orientation preferences de- 
rived from them are shown in Fig. [1^, b. Numerous studies confirmed that the orientation pref- 
erence of columns is an almost everywhere continuous function of their position in the cortex. 
The domains formed by neighboring columns with similar orientation preference are called iso- 
orientation domains 1251. 



1.2 Pinwheels 

At many locations the iso-orientation domains are arranged radially around a common center 
IISIlOl. Around these pin wheel ||TT| centers, stimulus orientations are represented in circular fash- 
ion. Such an arrangement had been previously hypothesized on the basis of electrophysiological 
experiments 1261 i27| and theoretical considerations [28 J. The regions exhibiting this kind of radial 
arrangement were termed pinwheels [11 J. The centers of pinwheels are point discontinuities of the 
field i9(x) where the mean orientation preference of nearby columns changes abruptly. They can 
be characterized by a topological charge 



which indicates in particular whether the orientation preference increases clockwise around the 
center of the pinwheel or counterclockwise. Here, Cj is a closed curve around a single pin wheel 
center at Xj. Since t9 is a cyclic variable within the interval [0,7r) and up to isolated points is a 
smooth function of x, g,; can only have the values 



where n is an integer number |[29ll . If its absolute value \qi\ is 1/2, each orientation is represented 
once in the vicinity of a pinwheel center. Pinwheel centers with a topological charge of ±1/2 
are simple zeros of 2;(x). In experiments only pinwheels that had the lowest possible topological 
charge qi = ±1/2 are observed. This means there are only two types of pinwheels: those whose 
orientation preference increases clockwise and those whose orientation preference increases coun- 
terclockwise. This organization has been confirmed in a large number of species and is therefore 
believed to be a general feature of visual cortical orientation maps |l30ll3Tll32ll33ll23ll34| . 

1.3 Hypercolumn and pinwheel density 

The pattern of preferred orientations is roughly repetitive |[7l,'2T]. The column spacing A, i.e. the 
spacing between adjacent iso-orientation domains preferring the same stimulus orientation, is typ- 
ically in the range of ~ 1mm. The column spacing A determines the size of the cortical hypercol- 
umn, which is considered to be the basic processing unit of the visual cortex IZll22l|355. We define 
the size of a hyper column by A^. The pinwheel density is defined as the number of pinwheels 
per unit area A^. Thus, by this definition, the pinwheel density is independent of the spacing of 
columns and dimension-less. 

1.4 Intra-cortical connectivity 

Visual cortical neurons are embedded in densely connected networks p6|. Besides a strong con- 
nectivity vertical to the cortical sheet between neurons from different layers within a column, neu- 
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rons also form extensive connections horizontal to the cortical surface linking different orientation 
columns. These connections extend for several millimeters parallel to the cortical surface and are 
therefore called long-range horizontal connections. As shown in Fig. [l}? for the example of the tree 
shrew, these connections are clustered primarily connecting domains of similar orientation prefer- 
ence. They have been observed in various mammals Il37l[23l l38l [121 and repeatedly hypothesized 
to play an important to role in visual processing tasks such as contour integration. 

1.5 Activity-dependent development 

In normal development, orientation columns first form at about the time of eye opening | l39l[l3l[T2l 
which for the ferret is approximately at post natal day (PD) 31. This is just a few days after neurons 
first respond to visual stimuli. A subset of neurons show orientation preference from that time on, 
but the adult pattern is not attained until seven weeks after birth |40 1. Roughly clustered horizontal 
connections are present by around PD 27 141] . Many lines of evidence suggest that the formation 
of orientation colunms is a dynamical process dependent on neuronal activity and sensitive to 
visual experience flTl |42|. This is suggested not only by the time line of normal development, 
but also receives support from various experiments manipulating the sensory input to the cortex. 
Most intriguingly, when visual inputs are rewired to drive what would normally become primary 
auditory cortex, orientation selective neurons and a pattern of orientation columns even forms in 
this brain region that would normally not at all be involved in the processing of visual information 
flSl l42||. This observation suggests that the capability to form a system of orientation columns 
is intrinsic to the learning dynamics of the cerebral cortex given appropriate inputs. Moreover, 
the comparison of development under conditions of modified visual experience demonstrates that 
adequate visual experience is essential for the complete maturation of orientation columns and that 
impaired visual experience, as with experimentally closed eye-lids [13, 12J or by rearing kittens in 
a striped environment consisting of a single orientation |,14J , can suppress or impair the formation 
of orientation columns (but see also |[43[| '). 

Viewed from a dynamical systems perspective, the activity-dependent remodeling of the cor- 
tical network is a process of dynamical pattern formation. In this picture, spontaneous symmetry 
breaking in the developmental dynamics of the cortical network underlies the emergence of cor- 
tical selectivities such as orientation preference |[T6[| . The subsequent convergence of the cortical 
circuitry towards a mature pattern of selectivities can be viewed as the development towards an 
attractor of the developmental dynamics flE]. This is consistent with the interpretation of cortical 
development as an optimization process. In the following, we will briefly describe a model |[19]| 
that is based on this view. 

1.6 Modeling cortical self-organization 

Self-organization has been observed to robustly produce large scale structures in various complex 
systems. Often, the class of patterns emerging depends on fundamental system properties such 
as symmetries rather than on system specific details. Pattern formation can therefore often be 
described by abstract models incorporating these properties only. For systems undergoing a Turing 
type instability, canonical model equations are of the S wif t-Hohenberg |[20l [2| type 

dtz = F[z] 

= LsHZ + N2[z] + N3[z] + --- (5) 



5 



where the linear part is 



LsH = r-(kl + VY (6) 



and z{-K,t) is a complex scalar field. If the bifurcation parameter r < 0, the homogeneous state 
2:(x) = is stable. For r > 0, a pattern with wavelength close to A = 2iT/kc emerges. The lowest 
order nonlinearities ^^^^2 and are quadratic and cubic in z, respectively. The form of these non- 
linearities determines the class of the emerging pattern, i.e. whether hexagons, rotating stripes, 
spiral waves, or another type of pattern emerges. 

Following this paradigm, we adopted a model of the form Eq. (|5]) with nonlinearities derived 
from key features of the visual cortex (following |[T9ll44[ ). As in experimental recordings, orienta- 
tion columns are represented by a complex field 1281 



z(x) = |z(x)|e^2^('') (7) 

where 'd is the orientation preference and \z\ a measure of selectivity at location x = (x, y) in the 
map. The factor 2 in the exponent accounts for the vr-periodicity of stimulus parameter orienta- 
tion. Constructing the nonlinearity of the model relies on the following assumptions. The model 
includes the effects of long-range intracortical connections between columns with similar orienta- 
tion preference (Fig. [l]:). Unlike many non-living systems in which interactions are purely local, 
long-range interactions are an important and distinctive feature of neuronal circuits in the cortex 
and particularly in the primary visual cortex. Further, based on the spatial homogeneity of circuits 
across cortex, it is assumed that the dynamics is symmetric with respect to translations, 

F[fy z] = fy F[z] with fy z(x) = z(x + y) , (8) 

and rotations 



F[Rf^ z] = Rp F[z] with R/j z(x) 



cos(/3) sin(/3) 
— sin(/3) cos(/?) 



(9) 



of the cortical sheet. This means that patterns that can be converted into one another by translation 
or rotation of the cortical layers belong to equivalent solutions of the model, Eq. (|5]) , by construc- 
tion. It is further assumed that the dynamics is symmetric with respect to shifts in orientation 

F[e''t' z] = e^'t' F[z] . (10) 

Thus, two patterns are also equivalent solutions of the model, if their layout of orientation domains 
is identical, but the preferred orientations differ everywhere by the same constant angle. Solutions 
shall contain representations of all stimulus orientations. For simplicity, couplings to other visual 
cortical representations such as ocular dominance or retinotopy are neglected. Considering only 
leading order terms up to cubic nonlinearities a nonlinearity fulfilling these requirements is given 
by 

N2[z{^)] = 

iV3[z(x)] = (l-5)k(x)|2z(x) + 

(9-2)1 dVi^<x(y-x)(z(x)|z(y)|2 + iz(x)z(y)2^ , (11) 
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with long-range interactions mediated through convolutions with a Gaussian 

KA^) = ^e-^ (12) 

with range a. The second parameter < g < 2 controls local and nonlocal interactions. The 
first term is the only strictly local term consistent with the required symmetries, the second non- 
local term represents the simplest non-local term that is symmetric with respect to these and with 
respect to permutations 



Ns{u,v,w) = Ns{w,u,v) . (13) 

Here the non-linear operator is written in a trilinear form as introduced in ||T9l l44ll . This addi- 
tional symmetry implies that all two-orientation solutions, for instance real valued solutions, of 
the model, are unstable, which in turn guarantees that all stimulus orientations are represented. 
One should note that in this model the spatial range of the nonlinearity a is a control parameter 
independent of the wavelength A. The patterns selected for different ratios a/ A are displayed in 
Fig. |2b. It can also be derived as an approximation to a model for the combined development of 
orientation preference and long-range horizontal connections BH . The model as defined in Eq. ^ 
[TSt is variational. It is consistent with synaptic models based on Hebbian plasticity, e.g. ||T7l[45ll46l . 
It is the only model known to the authors that exhibits stable aperiodic solutions dominated by a 
single spatial scale A. These solutions resemble orientation maps observed in the visual cortex. In 
the following, we discuss the structure of these solutions and the phase diagram of the model. 



1.7 Weakly nonlinear analysis 

Solutions of the model close to the bifurcation point r = are known in closed form, derived by 
means of a perturbation method called weakly nonlinear stability analysisflZlEl. When the dy- 
namics is close to a finite wavelength instability, the essential Fourier components of the emerging 
pattern are located on the critical circle. Solutions are planform patterns 

z(x) = J^A.e^'^^- (14) 

j 

composed of a finite number of Fourier components with wavevectors on the critical circle, | kj | = 
kc- By symmetry, the dynamics of amplitudes Ai of a planform are governed by amplitude equa- 
tions 

Ai = Ai — ^ ^ gij \Aj\ Ai — ^ ^ fijAjAj-A^- (15) 
j j 

where j~ denotes the index of the mode antiparallel to mode j. The form of Eq. dlSt is universal for 
models of a complex field z satisfying symmetry assumptions ((8l lT0]l . All model dependencies are 
included in the coupling coefficients gij and fij and may be obtained from F[z] by multiscale ex- 
pansion ||47l[2|. Denoting the angle between the wave vectors kj and kj by a and 6ij the Kronecker 
delta, the coefficients readlUlllll 

9ij = (^1 - ^^ij^ 9 (a) 

fij = {l-6i,-6,^^)f{a) (16) 
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Figure 2: a. Essentially complex planforms with different numbers n = 1, 2, 3, 5, 15 of active modes: 
The patterns of orientation preferences ^(x) are shown. The diagrams to the left of each pattern 
display the position of the wavevectors of active modes on the critical circle. For n = 3, there 
are two patterns; for n = 5, there are four; and for n = 15, there are 612 different patterns, b. 
Phase diagram of model. If non-local interactions are dominant {g < 1) and long-ranging (a large 
compared to A), quasiperiodic planforms are selected. Reproduced from |44| . 



where 

g{a) = g + 2(2 — g) exp (— cj^/c^) cosh (cr^/c^ cos(a)) 

/(«) = (17) 

are called angle-dependent interaction functions. 

Stationary solutions of Eq. flSl) are given by families of planforms lfT9l |44| 

n-l 

z(x) = ^|ylj|e^(^^-^^-^+'^^-^ (18) 

j=0 

of order n with wavevectors 




distributed equidistantly on the upper half of the critical circle and binary values Ij = ±1 deter- 
mining whether the mode with wave vector kj or with wavevector —hj is active. These planforms 
cannot realize a real valued function and are called essentially complex planforms (Fig. |2]|. For 
these planforms the third term in Eq. flSl) vanishes and the effective amplitude equations for the 
active modes reduce to a system of Landau equations 

A^ = Ai-Y,9ij\Ajf A (20) 
j 
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with stationary solutions Eq. ((181 ) with amplitudes of equal modulus 



) 



1/2 



A, 



(21) 



j 



and an arbitrary phase (pi independent of the mode configuration Ij . If the dynamics is stabilized 
by long-range nonlocal interactions {g < 1, a > A), large n planforms are the only stable solutions. 
In this long-range regime, the order of planforms grows as 



approximately linear with the interaction range. For a given order n, different planforms are de- 
generated in energy. This is a consequence of the permutation symmetry Eq. ([131 ). This symmetry 
also implies that the relevant stable solutions are essentially complex planforms which in turn 
guarantees that all stimulus orientations are represented. 

2 Calculation of pinwheel density 

2.1 Large range limit of interactions: Planform anisotropy 

First, we calculate the average pinwheel density pi for an ensemble of planforms, Eq. (|l8l) , with a 
fixed set of wavevector directions / = {Iq, h,. . . , In-i) but arbitrary phases 4>j in the limit n ^ oo. 
Here, and until noted otherwise, () shall denote average over phases In this limit, z and local 
linear functionals of z have Gaussian statistics such that the density of pinwheels is determined by 
the second order statistics of the field. Second, we evaluate the expectation value of pi over all sets 



For large a, to good approximation ga 1 and gij g and hence \Ai\ 1/ ^yng. Planforms ((18)> 
simplify to 



where for later convenience the constant ^/2g was absorbed into z(x). Pinwheels are the zeros of 
the field ^(x). The number of pinwheels in a given area A is obtained by 



n ~ 27ro"/A 



(22) 



oil. 




(23) 




(24) 



where 6{x) denotes Dirac's delta function and 



J(z(x)) 



ai?(x) gj(x) ai^(x) aj(x) 

dx dy dy dx 



(25) 



is the Jacobian of the field 

z(x) = i?(x) + i/(x) (26) 



split here for later convenience into its real and imaginary part. Averaging Eq. ((24l) over the 
ensemble of phases (pj reads 




9 



implying that 



PI = (5(z(x))J(z(x))) (28) 

is the expectation value of the pinwheel density for a fixed set of /. This expectation value only 
depends on local quantities, namely on the field, Eq. ([23l) , and its spatial derivatives 

Vz(x) = iJ-Yl ^ikje^('^'^^^+'^^) (29) 
^ j=o 

such that knowledge of the joint probability density p {z, Vz) is sufficient to evaluate Eq. ([28l) . 
Owing to the central limit theorem, this probability density becomes Gaussian in the large n limit. 
Eq. (|28l > is then determined by the first and second order statistics of z and Vz. Furthermore, since 
their statistics is the same at each location x, it is sufficient to evaluate Eq. ((28|) for z{0), Vz{0). The 
spatial dependency is thus omitted in the following. The average in Eq. ([28l) is given by an integral 
over the joint probability density 

p(v) = e-^'^^"''' (30) 

of components 

V = {R, I, d^R, dj, dyR, dyl) . (31) 

with covariance matrix C which shall be analyzed in the following. 

First, the diagonal elements of C are evaluated. Using dlSt and ([26l >, the auto-correlations of the 
field are 

2 

j,j'=o 

= 1 (32) 



for the real part and 



2 



) = -Yl (sm(/)j sm</)jv> 
jJ'=o 

= 1 (33) 
for the imaginary part. For the spatial derivatives one obtains 



2 



{dxR)l = — I jl ji kxjkxji (^sin sin 
^ j,j'=o 

^ n—l 

= (34) 



3=0 
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where k^j is the x-component of kj and, likewise, 

i 

!Wf) = ^E^^.- (35) 

i 

The equality of these correlations follows from inserting Eq. ( [191) into (|34]) yielding 



2 



j=0 



sm 

n 

j=0 



n 



2 n-l 



j=0 

\2 









\ n 




( jvr 




n 




( jtt' 




V n 



idyRY) (36) 



and also 



/ n ^ x"^ J 



j=0 

_ n-l 

i=o 

u2 //Q dn2 



2 



fc,^- (a,i?)^ (37) 



such that all auto-correlations become 



{d.Rf) = {{dyRf) = ({djf) = {{dylf) = f = 27r2 (38) 



when choosing without loss of generality the column spacing to be A = 27r/A;c = 1. 

Most off -diagonal elements of the covariance matrix C vanish. All non-vanishing contributions 
are related to the planform anisotropy defined by 



^ n— 1 



3=0 
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which depends on the set /. The covariance between the field and its derivative reads 



{zVz) 



j,j'=o 

n-l 



n 



j=0 



(40) 



with X = Ylj=o ^j^j — 0- The modulus X = Ixl is small for an isotropic distribution of wavevec- 
tors Ijkj. To estimate its upper bound Xmax, consider the most anisotropic case with all Ijhj situ- 
ated in the right plane {Ij = 1 for j < n/2, Ij = —1 for j > n/2). For large n this upper bound 
is 



Xn 



n-l 



3=0 



vr 



vr 



E 

7 = -- 
J 2 

2 

TT 

" 2 



TT 



n 



dae^ 



(41) 



such that the modulus ^ 



of the anisotropy is bounded by < ^ < 1. In the following, without 



loss of generality, we assume ^ = ^(1,0) implying that all correlations involving one derivative in 
y-direction vanish. Correlations involving dx are obtained by writing Eq. jlOl ) and 



n-l 



{zVz) 



n 



j,j'=o 







in the form 

{zVxz) = {RdxR) + {Idxl) + i {{IdxR) - {Rdxl)) = -iS^ 
(zVxz) = (Rd^R) - (IdJ) + i {{IdxR) + {RdJ)) = 

and comparing both imaginary parts showing that 

- {Id,R) = {RdJ) = 4C 



(42) 



(43) 



(44) 



does not vanish for anisotropic planforms. 

Expression ((44)l are the only non-vanishing non-diagonal elements of the matrix C. Indeed, 



{RdxR) = {IdJ) = 0, 



(45) 
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follows from comparing both real parts in Eq. ( [43l) . Furthermore, correlations between the real and 
imaginary part and between their derivatives, e.g. 



{RI) = 

(d^Rdj) = 
{dyRdj) = 



(46) 



vanish since they contain terms of the form(sin cpj cos (pj/ ) = (sin (pj ) (cos =0. Finally, because 

2 

(dxRdyR) = — Ijlj' kjxkj'y (sin (pj sin ) 

jj' 



> sm ivr— 

71 ^—^ \ n 



j=0 



C 



( {R^) {d.IR) 

(/') {dxRI) 

{IdxR) {{dxR?) 
{RdJ) {{djf) 



\ 

( 1 46 

1 -4?; 
-46 27r2 
46 27r2 



{{dyRf 



{{dyl?) J 



\ 



The integral we want to solve is 



27r2 / 



Pi 



(27r)3\/d^ 



For the exponent 



we find 



^ = — v^C" V 
2 



(47) 



vanishes for arbitrary n, also correlations between derivatives in different directions, 

(dxRdyR) = {djdyl) = , (48) 

do not contribute to the density of pinwheels. 

Altogether, the covariance matrix for the vector v = {R, /, dxR, dxl, dyR, dyl) reads 



(49) 



(50) 



(51) 
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Performing the integral over R and /, we get 



Pi 



1 / 

VdetC J 



d w J exp — - 
4 



7r2-8C2 



(27r)3Vdet 

where w = {dxR, dxl, dyR, dyl). Substituting 



+ 



where 



we have 

and the Jacobian reads 
such that we obtain 

Pi 



dxR = r2 008^2 

dyR = ricos^i 

dxl = r2sin02 

dyl = ri sin 9i 



< 01,^2 <27r 
< ri , r2 < oo 

d'^w = rir2dridr2d9id92 
J = rir2\ sm{di — 02)1 



dridr2rfr^ exp --( ^, + ^) 



(27r)2\/d^, 

after integrating over angles. Evaluating the determinant as 

det C = 167r^(7r2 - 8^2)2 
and performing the final integration steps, we obtain the pinwheel density 



Pi 



(52) 



(53) 



(54) 
(55) 

(56) 
(57) 

(58) 
(59) 



(60) 



of an ensemble of planforms with a fixed set of wavevector directions I. This result shows that 
the pinwheel density pi only depends on the anisotropy ^ of the considered planform. Since the 
anisotropy ranges within < ^ < 1, the pinwheel density is confined by 1.36 < < tt. 
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2.2 Distribution of planform anisotropics 

What is the distribution of anisotropies ^ in the large n limit? To address this question consider 
the ensemble defined by the different sets I = (/q, • • • ,^n-i)- From now on, () shall denote the 
expectation value over this ensemble. In the following we assume that the distribution of the 
vector anisotropies ^ is isotropic and Gaussian in the large n limit. The isotropy follows from the 
rotation symmetry of the model equations. It implies that 



0. 



(61) 



The assumption of Gaussian statistics of ^ is justified by the fact that by its definition ((39)) , the 
vector anisotropies ^ results from various wavevectors Zjkj with pairwise independent directions 



L. The distribution thus reads 



with variance given by 



P[0 = exp 



(62) 



IGn^ 



fi„ 2 '^J 



2 

j 

2 



4n 



(63) 



where {ijlj') = Sjji. The probability density for ^ follows from the distribution of the vector 
anisotropy ^ by 



P{0 = 27r^p 

8n f in . 



(64) 



where in the first equation the prefactor accounts for the change to polar coordinates. 
2.3 Pinwheel density in the large n limit 

Distribution (|64l ) and Eq. ( [60l) yield a distribution of pinwheel density pi by means of the coordinate 
transform,^ pi. With 



Pipi 



PiHPi)) 



dpi 



(65) 
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and ^ = v/ (vr^ — p^)/8 (Eq. ( |60l) ) the pinwheel density distribution reads 



P^pi") = ':;^P^ ~ ft ) j ' < ft < ^ • (66) 

The expectation value of the pinwheel density follows as 

/•TT 

{p) = / dpipip{pi) 







= TT ^^AV^ (67) 

where 5>j is the imaginary error function. In the limit n — > oo, the second term vanishes implying 

lim (p) = vr . (68) 

Moreover, the pinwheel density is (^-distributed at (p) = vr in the large n limit as its variance u 
converges towards 0. Having 

(p2> = r dpp^pip) 

Jo 

n 

2 2e 2 + n-2 

= TT'' (69) 

n 



which goes to tt^ for n ^ oo, one finds 



ip') - {Pf 

(70) 



suggesting that for large n not only the average pinwheel densities {p) but the density p of almost 
every realization is close to vr. Thus, in the limit of infinite interaction range, the value of the 
average pinwheel density is vr in almost every realization. 



2.4 Intermediate range of interaction 

Next, we ask to which extent these results do remain valid for finite interaction ranges and weather 
they are fairly robust against variations in range. Real cortices are of course finite and the effective 
interaction range for a given species or animal may depend on the range of intracortical long-range 
horizontal connections (see Fig. [l)c. It is therefore important to analyze the pinwheel density in 
solutions for intermediate interaction range a, i.e. in planforms of intermediate order n. To this 
end, we numerically synthesize planforms of various order n with randomly chosen Ij and (pj 
and calculate the pinwheel density averaged across planforms. Pinwheel centers are identified by 
the crossings of the zero contour lines of the real and imaginary part of the field z. To calculate 
the expectation value of the pinwheel density (p) in an ensemble of planforms of order n, we 
synthesize planforms (|T8]l with randomly chosen sets of wavevector directions Ij and phases (pj 
in a quadratic region with linear extension L on a 2048 x 2048 grid choosing an aspect ratio of 
F = L/A = 32. In each planform, pinwheel locations are identified within a quadratic subregion 
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Figure 3: Pinwheel densities for finite 
interaction range n « lira/ A. a. Aver- 
age pinwheel density (p) in planforms 
of order n (diamonds). Planforms with 
random sets of wavevector directions 
Ij and phases (pj were synthesized on 
2048 X 2048 grids with aspect ratio T = 
32. For each planform the pinwheel 
density was determined in a quadratic 
r = 8 subregion. Planforms were 
drawn until a achieving SEM A = 0.03 
for n < 15 and A = 0.01 for n > 15. The 
dashed line represents (p) = n, valid 
in the limit n — > 00. Note that the av- 
erage pinwheel density (p) is close to 
vr even at intermediate orders of n. h, 
Standard deviation (SD) s of densities 
p. Note that s decreases with order n. 



of size 8A x 8A. Realizations are collected until sufficient precision of average pinwheel densities 
(p) is reached, measured by the standard error measure (SEM) A = s/ \/]V, were s is the standard 
deviation (SD) of densities p and N the number realizations. The afforded precision of A < 0.03 
for 5 < n < 14 and A < 0.01 for 15 < n < 20 requires between 100 and 1000 realizations per 
order n. At a given order n, we determine the average pinwheel density (p) and the SD s from the 
ensemble of realizations. Note that by this method, both quantities can be evaluated with arbitrary 
precision by using sufficient large number of realizations. 

Fig. IHK shows the average pinwheel density (p) of planforms of various orders ranging between 
3 < n < 25. Among smaller orders n < 7, average densities (p) fluctuate substantially covering 
the range 2.5 < (p) < 3.5. However, they are much more confined, between 2.9 < (p) < 3.2, 
for intermediate orders 8 < n < 15. For large order n > 15, the ensemble average (p) appears 
to converge towards vr from below. This is consistent with Eq. (|60l ) showing that already for 
intermediate orders of n averaging over I leads to pinwheel densities smaller but close to vr. That 
average pinwheel densities (p) are smaller than vr is explained by the upper bound of pi in Eq. (|60l ) 
implying an upper bound also for its average. That they are, in fact, not much smaller than vr is 
suggested by the "relativistic" form of Eq. ( [6OI) . Even moderate anisotropies up to, e.g., £,1 = 0.3 
result in pinwheel densities (pi) > 3. Furthermore, as shown in Fig. [Sb, the variation s of pinwheel 
densities p in different realizations decreases drastically with n. For large n, s became successively 
smaller consistent with the limiting value = p = derived above. These results show that even 
for intermediate interaction ranges, both the average pinwheel density (p) of an ensemble of maps 
and the pinwheel density p of almost every single realization are close to vr. 
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3 Discussion 



We derived signatures of cortical self-organization that can be tested experimentally. In a model for 
the self-organization of the system of orientation preferences in the visual cortex we calculated the 
density of pinwheels (topological defects) that the model is predicting. We find that pinwheel den- 
sities close to vr are robustly selected if interactions between remote contour detectors are prevalent. 
Near criticality (r <^ 1), in the limit of large intracortical interaction range a, the average pinwheel 
density converges to the fixed number (p) = vr. For intermediate ranges lA < o" < 4 A, average 
pinwheel densities (p) remain smaller but close to this limit value. For successively larger ranges 
a, the pinwheel density approaches vr from below. Moreover, the variation of pinwheel densities p 
across realizations decreases with interaction range such that when the interaction range is large, 
almost every map exhibits a pinwheel density close to vr. Thus, for a broad parameter regime of 
interaction ranges, the model predicts an average pinwheel density close to vr even in individual 
cortical orientation maps. 

The results presented here are obtained close to the bifurcation point (r <^ 1). In this regime, 
the model is well approximated by the amplitude equations and its behavior is representative for 
a large class of models. Thus, the results derived in this regime are robust against variation of 
the model. Investigations further away from criticality at r > would be instructive, but can 
at present only be carried out by numerical methods. While dependencies on the details of the 
model are expected to play a role for r > 0, solutions in this regime may still be categorized by 
those obtained for r <^ I. 

In the visual cortex, distant neurons can directly and strongly interact by long-range horizon- 
tal connections and by feedback connections from one visual area to another. In mature cortical 
circuits, such connections link neurons that share similar receptive field properties and represent 
somewhat displaced locations in visual space Il37ll23l . In fact, in a wide range of mammals, the 
range of these connections is much larger than the spacing of columns 1371 l23l [12J. Thus, these 
species are candidates for testing the value of the average pinwheel density predicted by our study. 

The approach presented in this study, to our knowledge, for the first time allows to analytically 
characterized the statistics of pinwheels in populations of asymptotically stable states of visual 
cortical self-organization. The results obtained are valid for the general class of Turing-type sys- 
tems exhibiting the four symmetries Eq. ((8] |10I13|| and dominated by long-ranging interactions 
when studied close to the instability threshold. Asymptotically stable states have been studied in 
two competing model classes. In one class of models, pinwheels annihilate dynamically, such that 
pinwheel density is a time dependent quantity and no particular finite pinwheel density is intrin- 
sically selected IIT8ll48ll49ll50| . In other models, a substantial number of pinwheels are preserved in 
the final state, but the pinwheels typically crystallize into repetitive arrangements l|48ll5Tl[49ll52l . 
These states are easy to characterize, but are not even qualitatively consistent with the aperiodic 
arrangement of visual cortical orientation columns. The model studied here is distinguished from 
these competing models by exhibiting in the long-range interaction dominated regime {g < 1, 
a > A) a large multiplicity of aperiodic stable solutions. Our study presents a transparent an- 
alytical approach for studying pinwheel density selection in such models. It will be interesting 
to generalize this approach to further model classes to comprehensively clarify whether genuinely 
different models of visual cortical development can be distinguished through their pinwheel statis- 
tics, as the above results suggest. 

It is natural to ask whether and how the predicted pinwheel density can be tested experimen- 
tally. The currently available published data is consistent with a pinwheel density around 3 but 
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appears insufficient to obtain a high precision estimate of this quantity (see e.g. 1531 |54l, for a 
summary of the earlier literature see |18|). To obtain a precise estimate both pattern wavelength 
and absolute pinwheel number need to be reliably quantified. Over the past years, image analy- 
sis methods have been devised that enable estimating the pattern wavelength of orientation maps 
with a precision in the range of a few percent ( |[35l [55]| ). For instance, using a wavelet method 
for the estimation of local pattern wavelength Kaschube and coworkers showed that genetically 
related cats often differ in the mean wavelength of orientation columns by less than 4 % (|35|). As 
the count variance of pinwheel number estimates is of order (number of pinwheels)^/^ , estimating 
the pinwheel density with an accuracy in the percent range will require data sets encompassing on 
the order of 10000 genuine pin wheels. In many model animals, this size of data set is equivalent 
to more 100 brain hemispheres that need to be imaged under consistent experimental conditions. 
Whereas this may appear a large animal cohort for any individual study, such large size data sets 
will eventually accumulate in many laboratories that use intrinsic signal optical imaging as a stan- 
dard technique, enabling to test models of visual cortical development with quantitative precision. 
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